New Normal - Shifting Climate Baselines
Adam A. Kemberling
Updated on: 2021-05-31
# Turn on thematic for GMRI styling of plots:
# Set Style options
# Turn on thematic:
# Set Style options
# thematic_on(
# bg = NA,
# fg = NA,
# accent = NA,
# font = "auto"
# )
# Turn on the gmri font for plots - doesn't really connect to gmri font
showtext::showtext_auto()Shifting Norms - Moving Climate Norms to the New Climatology Period
# load average global climatologies as annual averages
yr_avgs <- str_c(box_paths$oisst, "daily_climatologies/yearly_means/")
clim_82_avg <- stack(str_c(yr_avgs, "avg_clim_1982to2011.nc "))
clim_91_avg <- stack(str_c(yr_avgs, "avg_clim_1991to2020.nc"))
clim_shift_avg <- stack(str_c(yr_avgs, "avg_clim_shift_82to91baselines.nc"))
# climatology_lists
grid_list <- list(
"1982-2011" = clim_82_avg,
"1991-2020" = clim_91_avg,
"Climate Shift" = clim_shift_avg)
# rotate them
grid_list <- map(grid_list, raster::rotate)# convert annual_avgs to stars objects
clim_82_st <- st_as_stars(grid_list$`1982-2011`)
clim_91_st <- st_as_stars(grid_list$`1991-2020`)
clim_shift_st <- st_as_stars(grid_list$`Climate Shift`)
# Set the crs to world molleweide
world_moll <- st_crs("+proj=moll")
# Warp to world projection
clim_82_moll <- warp_grid_projections(clim_82_st, projection_crs = "world moll")
clim_91_moll <- warp_grid_projections(clim_91_st, projection_crs = "world moll")
clim_shift_moll <- warp_grid_projections(clim_shift_st, projection_crs = "world moll")
# Get temp limits that fit both
overall_max <- max(c( max(clim_82_st[[1]], na.rm = T) , max(clim_91_st[[1]], na.rm = T) ))
overall_min <- max(c( min(clim_82_st[[1]], na.rm = T) , min(clim_91_st[[1]], na.rm = T) ))
temp_limits <- c(overall_min, overall_max)
# Map the two periods
plot_period <- function(in_stars_obj, period_label, crs_projection){
period_plot <- ggplot() +
geom_stars(data = in_stars_obj) +
geom_sf(data = st_transform(world, crs_projection), fill = "gray80", size = 0.2) +
scale_fill_distiller(palette = "RdYlBu", na.value = "transparent", limit = temp_limits) +
guides("fill" = guide_colorbar(title = expression("Average Annual Sea Surface Temperature"~~degree~C),
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
coord_sf(crs = crs_projection) +
map_theme +
labs(title = period_label)
return(period_plot)
}# Plotting the two periods
plot_82 <- plot_period(clim_82_moll, period_label = "1982-2011 Climatology", crs_projection = world_moll)
plot_91 <- plot_period(clim_91_moll, period_label = "1991-2020 Climatology", crs_projection = world_moll)
# Plotting them stacked
stacked_clims <- (plot_82 + theme(legend.position = "none")) / plot_91
stacked_clims
#### Shift in Climatology Plot
# Set diverging limits for scale
shift_limits <- max(abs(clim_shift_st[[1]]), na.rm = T) * c(-1,1)
# Plot the shift in the two
plot_shift <- ggplot() +
geom_stars(data = clim_shift_moll) +
geom_sf(data = st_transform(world, world_moll), fill = "gray80", size = 0.2) +
scale_fill_distiller(palette = "RdYlBu", na.value = "transparent", limit = shift_limits) +
guides("fill" = guide_colorbar(title = expression("Sea Surface Temperature Difference"~~degree~C),
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
map_theme +
coord_sf(crs = world_moll) +
labs(title = "Global Shifts in Annual SST from Climatology Transition")# Plotting the change
plot_shift
Region Extents
The following figures look specifically at two regions. 1.) The Gulf of Maine, and 2.) the Northeastern United States’ Continental Shelf. Both of these areas have been mapped below for clarity on what data is included for any regional metric.
# load shapes/timeseries paths
gmri_areas <- get_timeseries_paths(region_group = "gmri sst focal areas")
nelme_regions <- get_timeseries_paths(region_group = "nelme regions")
# get region shapes
gom <- read_sf(gmri_areas$apershing_gulf_of_maine$shape_path)
ne_shelf <- read_sf(nelme_regions$NELME$shape_path)
# lat/lon limits
map_lims <- list(x = c(-78, -63), y = c(35, 47))
# Gulf of Maine
gom_map <- base_map +
geom_sf(data = gom,
alpha = 0.4,
color = gmri_cols("gmri blue"),
fill = gmri_cols("gmri blue")) +
scale_y_continuous(breaks = seq(30, 50, by = 1)) +
scale_x_continuous(breaks = seq(-180, 180, by = 5)) +
coord_sf(xlim = map_lims$x, ylim = map_lims$y) +
map_theme +
labs(title = "Gulf of Maine")
#NE Shelf
shelf_map <- base_map +
geom_sf(data = ne_shelf, alpha = 0.4,
color = gmri_cols("teal"),
fill = gmri_cols("teal")) +
scale_y_continuous(breaks = seq(30, 50, by = 1)) +
scale_x_continuous(breaks = seq(-180, 180, by = 5)) +
coord_sf(xlim = map_lims$x, ylim = map_lims$y) +
map_theme +
labs(title = "Northeast Shelf")
# Mapping Them together
gom_map | shelf_map
Change in Annual Regional Temperatures
# masking function
mask_sf <- function(in_ras, sf_mask){
r1 <- crop(in_ras, sf_mask)
r2 <- mask(r1, sf_mask)
return(r2)}
# Mask them all
masked_grids <- map(grid_list, function(in_ras){
# Get means
gom_masked <- mask_sf(in_ras, gom)
shelf_masked <- mask_sf(in_ras, ne_shelf)
# return both as list
return(list(
"Gulf of Maine" = gom_masked,
"Northeast Shelf" = shelf_masked
))
})Gulf of Maine Climate Shift
# Pull Stars Object
gom_shift_st <- st_as_stars(masked_grids$`Climate Shift`$`Gulf of Maine`)
# Plot the shift in Climate
ggplot() +
geom_stars(data = gom_shift_st) +
geom_sf(data = new_england, fill = "gray80", size = 0.2) +
geom_sf(data = canada, fill = "gray80", size = 0.2) +
geom_sf(data = greenland, fill = "gray80", size = 0.2) +
scale_y_continuous(breaks = seq(30, 50, by = 1)) +
scale_x_continuous(breaks = seq(-180, 180, by = 5)) +
coord_sf(xlim = c(-71.25, -65), ylim = c(41, 44.5)) +
scale_fill_distiller(palette = "RdYlBu", na.value = "transparent", limit = shift_limits) +
guides("fill" = guide_colorbar(title = expression("Sea Surface Temperature Difference"~~degree~C),
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
map_theme +
labs(title = "Gulf of Maine Difference in Climate")
Northeast Shelf Climate Shift
# Pull Stars Object
ne_shift_st <- st_as_stars(masked_grids$`Climate Shift`$`Northeast Shelf`)
# Plot the shift in Climate
ggplot() +
geom_stars(data = ne_shift_st) +
geom_sf(data = new_england, fill = "gray80", size = 0.2) +
geom_sf(data = canada, fill = "gray80", size = 0.2) +
geom_sf(data = greenland, fill = "gray80", size = 0.2) +
scale_y_continuous(breaks = seq(30, 50, by = 1)) +
scale_x_continuous(breaks = seq(-180, 180, by = 5)) +
coord_sf(xlim = map_lims$x, ylim = map_lims$y) +
scale_fill_distiller(palette = "RdYlBu", na.value = "transparent", limit = shift_limits) +
guides("fill" = guide_colorbar(title = expression("Sea Surface Temperature Difference"~~degree~C),
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
map_theme +
labs(title = "Northeast Shelf Difference in Climate")
Annual Changes
# Get Average Temp and Change in Climate
shift_table <- imap_dfr(masked_grids, function(in_ras, climatology){
# Get means
gom_mean <- in_ras[["Gulf of Maine"]] %>% cellStats(mean, na.rm = T)
ne_mean <- in_ras[["Northeast Shelf"]] %>% cellStats(mean, na.rm = T)
# build table
grid_shift <- data.frame(
Climatology = rep(climatology, 2),
region = c("Gulf of Maine", "Northeast Shelf"),
avg_temp = c(gom_mean, ne_mean))
return(grid_shift) }) %>%
mutate(Climatology = factor(Climatology, levels = c("Climate Shift", "1991-2020", "1982-2011")))
# Plot the shifts
ggplot(shift_table, aes(x = avg_temp, y = Climatology)) +
geom_segment(aes(xend = 0, yend = Climatology), color = gmri_cols("gmri blue")) +
geom_point(size = 2, color = gmri_cols("orange")) +
facet_wrap(~region) +
geom_vline(xintercept = 0, color = "gray30", linetype = 3) +
labs(y = "", x = expression("Average Temperature"~~degree~C)) 
Changes to Seasonality
Loading Climatology and Heatwave Info
# Load regional timeseries
gom_ts <- read_csv(gmri_areas$apershing_gulf_of_maine$timeseries_path, col_types = cols()) %>%
mutate(time = as.Date(time))
ne_shelf_ts <- read_csv(nelme_regions$NELME$timeseries_path, col_types = cols()) %>%
mutate(time = as.Date(time))
# Side by side - Heatwave Detection
# Run climatology for both periods, merge and compare
hw_comparison <- function(region_timeseries){
# Run heatwave detection using new climate period
heatwaves_82 <- pull_heatwave_events(region_timeseries,
threshold = 90,
clim_ref_period = c("1982-01-01", "2011-12-31")) %>%
select(time, sst, seas_82 = seas, anom_82 = sst_anom, hw_thresh_82 = mhw_thresh, cs_thresh_82 = mcs_thresh,
hwe_82 = mhw_event, cse_82 = mcs_event, hwe_no_82 = mhw_event_no, cse_no_82 = mcs_event_no)
# Run heatwave detection using new climate period
heatwaves_91 <- pull_heatwave_events(region_timeseries,
threshold = 90,
clim_ref_period = c("1991-01-01", "2020-12-31")) %>%
select(time, sst, seas_91 = seas, anom_91 = sst_anom, hw_thresh_91 = mhw_thresh, cs_thresh_91 = mcs_thresh,
hwe_91 = mhw_event, cse_91 = mcs_event, hwe_no_91 = mhw_event_no, cse_no_91 = mcs_event_no)
# Subtract old from the new - set up flat date for plots
base_date <- as.Date("2000-01-01")
heatwaves_out <- left_join(heatwaves_82, heatwaves_91, by = c("time", "sst")) %>%
mutate(anom_shift = anom_91 - anom_82,
clim_shift = seas_91 - seas_82,
upper_shift = hw_thresh_91 - hw_thresh_82,
lower_shift = cs_thresh_91 - cs_thresh_82,
hwe_change = case_when(
hwe_82 == T & hwe_91 == F ~ "No Longer HW",
hwe_82 == F & hwe_91 == T ~ "New HW",
hwe_82 == T & hwe_91 == T ~ "Still HW",
hwe_82 == F & hwe_91 == F ~ "Still No HW",
T ~ "missing case"),
shift_direction = ifelse(clim_shift >= 0, "Temperature Increase", "Temperature Decrease"),
year = year(time),
yday = yday(time),
flat_date = as.Date(yday-1, origin = base_date))
return(heatwaves_out)
}
# Run comparisons
gom_comp <- hw_comparison(region_timeseries = gom_ts)
ne_comp <- hw_comparison(region_timeseries = ne_shelf_ts)Plotting Seasonality
# Plotting the seasonal variations
seasonality_check <- function(comp_data){
# Pull climatology
seasonal_pattern <- comp_data %>%
filter(between(time, as.Date("2019-01-01"), as.Date("2019-12-31"))) %>%
distinct(flat_date, .keep_all = T)
# Plot 82 clim
p_82 <-ggplot(seasonal_pattern, aes(x = time)) +
geom_line(aes(y = seas_82),
color = gmri_cols("gmri blue"),
size = 1) +
geom_line(aes(y = hw_thresh_82), linetype = 2, color = "gray60") +
geom_line(aes(y = cs_thresh_82), linetype = 2, color = "gray60") +
labs(x = "", y = expression("SST"~~degree~C),
subtitle = "1982-2011 Climatology") +
scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0))
# Plot 91 Clim
p_91 <- ggplot(seasonal_pattern, aes(x = time)) +
geom_line(aes(y = seas_91),
color = gmri_cols("teal"),
size = 1) +
geom_line(aes(y = hw_thresh_91), linetype = 2, color = "gray60") +
geom_line(aes(y = cs_thresh_91), linetype = 2, color = "gray60") +
labs(x = "", y = expression("SST"~~degree~C),
subtitle = "1991-2020 Climatology") +
scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0))
# # Plot both lines together
# p_both <- ggplot(seasonal_pattern, aes(x = time)) +
# geom_ribbon(aes(ymin = seas_82, ymax = seas_91), fill = "gray80") +
# geom_line(aes(y = seas_82),
# color = "1982-2011 Climatology",
# size = .75) +
# geom_line(aes(y = seas_91),
# color = "1991-2020 Climatology",
# size = .75) +
# scale_colour_manual(
# values = c("1982-2011 Climatology" = as.character(gmri_cols("gmri blue")),
# "1991-2020 Climatology" = as.character(gmri_cols("teal")))) +
# labs(x = "", y = expression("SST"~~degree~C),
# subtitle = "Expected Temperature from Climatologies") +
# scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0))
# Plot the difference
p_diff <- ggplot(seasonal_pattern, aes(x = time, y = clim_shift)) +
geom_segment(aes(yend = 0, xend = time, color = shift_direction)) +
scale_color_gmri(reverse = T) +
labs(x = "",
y = expression("Temperature Shift"~~degree~C),
color = "Shift in Climate Norm",
subtitle = "Difference in Climatologies") +
scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0)) +
guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) +
theme(legend.position = "bottom")
return(
list(
"clim82" = p_82,
"clim91" = p_91,
#"clim_both" = p_both,
"clim_shift" = p_diff)
)
}Gulf of Maine
# Return the plots
gom_seasonality <- seasonality_check(gom_comp)
# Assemble stack
(gom_seasonality[[1]] + labs(title = "Gulf of Maine")) / gom_seasonality[[2]] / gom_seasonality[[3]]
Northeast Shelf
# Return the plots
ne_seasonality <- seasonality_check(ne_comp)
# Assemble stack
(ne_seasonality[[1]] + labs(title = "Northeast Shelf")) / ne_seasonality[[2]] / ne_seasonality[[3]]
Heatwave Events
Heatmap of Marine Heatwave Events
# Plotting how anomalies are now suppressed
heatwave_checks <- function(comp_data){
# Set palette limits to center it on 0 with scale_fill_distiller
limit <- max(abs(comp_data$anom_91)) * c(-1,1)
# Base date for rectangle fill for NA
base_date <- as.Date("2000-01-01")
# Assemble heatmap plot
heatwave_heatmap <- ggplot(comp_data, aes(x = flat_date, y = year)) +
# background box fill for missing dates
geom_rect(xmin = base_date,
xmax = base_date + 365,
ymin = min(comp_data$year) - .5,
ymax = max(comp_data$year) + .5,
fill = "gray30",
color = "transparent") +
# Tile for sst colors
geom_tile(aes(fill = anom_91)) +
# Points for heatwave events
geom_point(data = filter(comp_data, (hwe_82 == TRUE | hwe_91 == TRUE)),
aes(x = flat_date,
y = year,
color = hwe_change),
size = .25) +
scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0)) +
scale_y_continuous(limits = c(1980.5, 2021.5), expand = c(0,0)) +
# Colors
scale_fill_distiller(palette = "RdYlBu", limit = limit) + # for actual anomalies
scale_color_manual(values = c("gray70", "black")) +
#5 inches is default rmarkdown height for barheight
guides("fill" = guide_colorbar(title = "Anomaly from 91-2020 Climatology",
title.position = "top",
title.hjust = 0.5,
barwidth = unit(3, "inches"),
frame.colour = "black",
ticks.colour = "black"),
color = guide_legend(title = "Change in Heatwave Record",
title.position = "top",
title.hjust = 0.5,
override.aes = list(size = 3))) +
labs(x = "",
y = "") +
theme_classic() +
theme(legend.position = "bottom",
panel.border = element_rect(color = "black", fill = "transparent", size = 1))
# Assemble pieces
return(heatwave_heatmap)
}Gulf of Maine
heatwave_checks(gom_comp) +
labs(title = "Shifting Heatwave Baseline - Gulf of Maine")
Northeast Shelf
heatwave_checks(ne_comp) +
labs(title = "Shifting Heatwave Baseline - Northeast Shelf")
Impact on Heatwave Events
Characteristics for Each Climatology
#### Summarize Yearly Heatwave Totals and Characteristics
hw_characteristics <- function(comparison_data, clim_group){
# Switch between 81 and 92
event_number <- switch(clim_group,
"82" = sym("hwe_no_82"),
"91" = sym("hwe_no_91"))
event_flag <- switch(clim_group,
"82" = sym("hwe_82"),
"91" = sym("hwe_91"))
sst_anom <- switch(clim_group,
"82" = sym("anom_82"),
"91" = sym("anom_91"))
clim_period <- switch(clim_group,
"82" = "1982-2011",
"91" = "1991-2020")
# 1. Characterize each heatwave event:
# number of heatwaves
# average heatwave duration
# remove NA as a distinct heatwave number
wave_events <- comparison_data %>%
mutate(yr = year(time)) %>%
group_by(yr, !!event_number) %>%
summarise(total_days = sum(!!event_flag, na.rm = T),
avg_anom = mean(!!sst_anom, na.rm = T),
peak_anom = max(!!sst_anom, na.rm = T),
.groups = "keep") %>%
ungroup() %>%
drop_na()
# 2. Get Annual Means:
wave_summary <- wave_events %>%
group_by(yr) %>%
summarise(num_waves = n_distinct(!!event_number),
avg_length = mean(total_days, na.rm = T),
avg_intensity = mean(avg_anom, na.rm = T),
peak_intensity = max(peak_anom, na.rm = T),
.groups = "keep") %>%
ungroup() %>%
mutate(clim = clim_period)
return(wave_summary)
}
#### Annual Heatwave Summary Details
gom_hw_82 <- hw_characteristics(comparison_data = gom_comp, clim_group = "82")
gom_hw_91 <- hw_characteristics(comparison_data = gom_comp, clim_group = "91")
ne_hw_82 <- hw_characteristics(comparison_data = ne_comp, clim_group = "82")
ne_hw_91 <- hw_characteristics(comparison_data = ne_comp, clim_group = "91")Shift in Heatwave Characteristics
# Join the two periods together for climatology shift plotting, doesn't overwrite
join_hw_summs <- function(summ_82, summ_91){ # Put the data together
summ_91 <- select(summ_91,
yr,
w = num_waves,
x = avg_length,
y = avg_intensity,
z = peak_intensity)
# Need to join to preserve all years
both_details <- summ_82 %>%
select(-clim) %>%
full_join(summ_91, by = "yr") %>%
mutate(w = ifelse(is.na(w), 0, w),
x = ifelse(is.na(x), 0, x),
y = ifelse(is.na(y), 0, y),
z = ifelse(is.na(z), 0, z))
return(both_details)
}
# Run the summary for both areas
gom_shift <- join_hw_summs(gom_hw_82, gom_hw_91)
ne_shift <- join_hw_summs(ne_hw_82, ne_hw_91)Plotting Heatwave Characteristic Shifts
#### Plotting shifts using dumbbell plot
hw_dumbbell <- function(wave_summary){
# 1. number of heatwaves
hw_counts <- ggplot(wave_summary) +
geom_segment(aes(x = yr,
y = num_waves,
xend = yr,
yend = w),
color = "gray60",
lineend = "round",
linejoin = "mitre",
size = 0.5,
arrow = arrow(length = unit(0.25, "cm"))) +
geom_point(aes(x = yr,
y = num_waves,
color = "1982-2011 Climatology"), size = 2) +
geom_point(aes(x = yr,
y = w,
color = "1991-2020 Climatology"), size = 2) +
geom_hline(yintercept = 0, color = "gray30", linetype = 3) +
labs(y = "# of Events",
x = "",
color = "") +
scale_color_gmri()
# 2. Average duration
hw_lengths <- ggplot(wave_summary) +
geom_segment(aes(x = yr,
y = avg_length,
xend = yr,
yend = x),
color = "gray60",
lineend = "round",
linejoin = "mitre",
size = 0.5,
arrow = arrow(length = unit(0.25, "cm"))) +
geom_point(aes(y = avg_length,
x = yr,
color = "1982-2011 Climatology"), size = 2) +
geom_point(aes(x = yr,
y = x,
color = "1991-2020 Climatology"), size = 2) +
geom_hline(yintercept = 0, color = "gray30", linetype = 3) +
labs(y = "Event Length",
x = "",
color = "") +
scale_color_gmri()
# 3. avg temp
hw_temps <- ggplot(wave_summary) +
geom_segment(aes(x = yr,
y = avg_intensity,
xend = yr,
yend = y),
color = "gray60",
lineend = "round", # See available arrow types in example above
linejoin = "mitre",
size = 0.5,
arrow = arrow(length = unit(0.25, "cm"))) +
geom_point(aes(x = yr,
y = avg_intensity,
color = "1982-2011 Climatology"), size = 2) +
geom_point(aes(x = yr,
y = y,
color = "1991-2020 Climatology"), size = 2) +
geom_hline(yintercept = 0, color = "gray30", linetype = 3) +
labs(x = "",
y = "Avg. Anomaly",
color = "") +
scale_color_gmri()
# 4. peak temp
hw_peaks <- ggplot(wave_summary) +
geom_segment(aes(x = yr,
y = peak_intensity,
xend = yr,
yend = z),
color = "gray60",
lineend = "round",
linejoin = "mitre",
size = 0.5,
arrow = arrow(length = unit(0.25, "cm"))) +
geom_point(aes(x = yr,
y = peak_intensity,
color = "1982-2011 Climatology"), size = 2) +
geom_point(aes(x = yr,
y = z,
color = "1991-2020 Climatology"), size = 2) +
geom_hline(yintercept = 0, color = "gray30", linetype = 3) +
labs(x = "Year",
y = "Peak Anomaly",
color = "") +
scale_color_gmri()
# Arrange them
column_plots <- (hw_counts / hw_lengths / hw_temps / hw_peaks) + plot_layout(guides = "collect") & theme(legend.position = "bottom")
return(column_plots)
}Gulf of Maine
hw_dumbbell(gom_shift) + plot_annotation(title = "Gulf of Maine - Shifting Baselines")
Northeast Shelf
hw_dumbbell(ne_shift) + plot_annotation(title = "Northeast Shelf - Shifting Baselines")